###############################################################################
# sim.function.R
#
# wrapper function to be called by snow.
# contains coalitionsim.R, fuel.R, drive.R, makedatasets.R from Gijs original paper
#
# Last edited: Jan 26 2015 by Roni 
###############################################################################

load.make.path <- function(x, folder, experiment) {
    load(paste0(folder, "/sim_",experiment,"_", x, ".Rdata"))
    return(as.data.table(dataset))
}

wrapper.sim.function <- function(input.list, voters, folder, potential.rules) {
  experiment <- input.list$experiment[1]
  end <- nrow(input.list)
 if (file.exists(paste0(folder, "/sim_",experiment, ".rds"))==FALSE) {
  apply(input.list, 1, sim.function,
        	voters=voters,
	        folder=folder,
	        potential.rules=potential.rules)
	   dat <- lapply(1:end, load.make.path, folder=folder, experiment=experiment)
	   dat <- rbindlist(dat)
	   saveRDS(dat, paste0(folder, "/sim_",experiment, ".rds"))
	   lapply(1:end,
	                function(x) file.remove(paste0(folder, "/sim_",experiment,"_", x, ".Rdata")))
 }
}



sim.function <- function(input, voters, folder, potential.rules) {
 num <- as.numeric(input[11])
 experiment <- as.numeric(input[7])

if (file.exists(paste0(folder, "sim_",experiment,"_", num, ".Rdata"))==FALSE) {
    
#test if other CPU is already working on this simulation
 set.seed(as.integer(input[10])) #make replicatable
 parties <- as.numeric(input[1])
 discount <- as.numeric(input[2])
 memory <- as.numeric(input[3])
 aspiration <- as.numeric(input[4])
 alpha <- rep(as.numeric(input[5]), parties)
 sd.distance <- as.numeric(input[6])*.12
 no.evo.length <- as.numeric(input[8])
 evo.burn.in <- as.numeric(input[9])
 run.seed <- as.numeric(input[10])
 start.rules <- sample(potential.rules, parties, TRUE)
 
 ################
 # COALITIONSIM #
 ################

 #get values for probability choice in first rounds (or if rule not used)
  runs <- no.evo.length+evo.burn.in+100 
  #parties <- length(start.rules)
  speed <- matrix(0.012,1,parties) 
  #discount <- round(runif(1),2) #discount of caretaker government= old.util * (1-discount)
  #memory <- 2
  rule <- matrix(NA, runs, parties)
  nrules <- 10 ## even if not 10 rule in simulation, rule label (1,2,3..) remain
  rule[1:runs,] <- rep(start.rules, each=runs)
  
 #rule[11,1] <- 1; rule[11,2] <- 2; rule[11,3] <- 3 #new rules not included yet  ## New code below 
 #if(parties>3){rule[11,4:parties] <- sample(c(1:3), (parties-3), replace=TRUE)} #new rules not included yet ## New code below
 #for (p in 1:parties){rule[12:20,p] <- rule[11,p]} ##New code below

 
 ########
 # FUEL #
 ########
 # set voter preferences
 vprefdim1 <- rtnorm(voters, (sqrt(2)/4), .12, lower=0, upper=(sqrt(2)/2))
 vprefdim2 <- rtnorm(voters, (sqrt(2)/4), .12, lower=0, upper=(sqrt(2)/2))
 mean.vprefdim1 <- mean(vprefdim1)
 mean.vprefdim2 <- mean(vprefdim2)
 vpref <- cbind(vprefdim1, vprefdim2)
 gmedian <- l1median(vpref)
 sum.dist.voter.median.voter <- sum(apply(vpref, 1, euclidean.distancer,
                                          gmedian[1], gmedian[2]))


 ## #rescale voter preferences to the interval [0,10]
 ## vprefdim1 <- as.matrix(rescale(vprefdim1))
 ## vprefdim2 <- as.matrix(rescale(vprefdim2))

 distances <- matrix(NA, voters, parties)
 vote <- matrix(NA, voters, 1)
 representation <- matrix(NA, runs, 1)
 gov.representation <- matrix(NA, runs, 1)
 eccentricity <- matrix(NA, runs, parties)
 
 # set information storages
 elections <- matrix(NA, runs, parties)
 supdim1 <- matrix(NA, runs, parties)
 supdim2 <- matrix(NA, runs, parties)
 
 # set party position
 partydim1 <- matrix(NA, runs, parties)
 partydim2 <- matrix(NA, runs, parties)
 govdim1 <- matrix(NA, runs, 1)
 govdim2 <- matrix(NA, runs, 1)
 govseats <- matrix(NA, runs, parties)
 caretaker <- matrix(NA, runs, 1)

 pprefdim1 <- rtnorm(parties, (sqrt(2)/4), sd.distance, lower=0, upper=(sqrt(2)/2)) #parties' policy preferences
 pprefdim2 <- rtnorm(parties, (sqrt(2)/4), sd.distance, lower=0, upper=(sqrt(2)/2)) #parties' policy preferences

 partydim1[1,] <- pprefdim1
 partydim2[1,] <- pprefdim2

 # Information
 util <- matrix(NA, runs, parties)
 bargpow <- matrix(NA, runs, parties) 
 bestrule1 <- matrix(NA, runs, nrules)
 bestrule2 <- matrix(NA, runs)
 votpref1 <- matrix(0, runs, 2)
 votpref2 <- matrix(0, runs, 2)
 sysdens <- matrix(NA, runs, 1)
 enr <- rep(NaN, runs)
 vote.sums <- matrix(NaN, nrow=runs, ncol=10)



 distances  <- t(apply(cbind(vprefdim1, vprefdim2), 1, euclidean.distancer, partydim1[1,],
                partydim2[1,]))



 vote <- apply(distances, 1, voter) #proximity voting
   #how far are voters, on avarage, from their closest party away in this party system?

  elections[1,] <- tabulate(vote, parties)
 
 # Government formation
 gov <- matrix(NaN, runs, parties)
 valid <- 0
 while (valid==0) {#get a incumbent government with at least one party that has seats as member
  gov[1,] <- as.vector(rbinom(parties, 1, .5))
  valid <- ifelse(sum(gov[1,]*elections[1,])==0,0,1)
 }
 gov[2,] <- gov[1,]

 #I exclude the first government because it contains no parties.
 allc <- expand.grid(rep(list(c(0,1)),parties))[-1,] #get all possible governments
 rownames(allc) <- NULL #because first government is excluded rownames need to be reset
 colnames(allc) <- NULL #to allow for comparisons with identical


 ## get utility for parties given randomly chosen government and its position
  util[1,] <- initial.utility(vprefdim1, vprefdim2, gov[1,], partydim1[1,], partydim2[1,], alpha,
                             pprefdim1, pprefdim2)

  seatsc <- ifelse(gov[1,]==1, elections[1,], 0)
  sumseats <- sum(seatsc)
  govdim1[1,] <- sum(partydim1[1,]*seatsc/sumseats)
  govdim2[1,] <- sum(partydim2[1,]*seatsc/sumseats)
  govseats[1,] <- ifelse(gov[1,]==1, elections[1,], 0)

 

 p <- seq_len(parties)


 print(c("r=", 2))
  # Party positioning for evolution

 
   ## #get position of government party with highest bargpower  
   ## dim1.largest.gov <- partydim1[(r-1),][bargpow[(r-1),]==max(bargpow[(r-1),])]
   ## dim2.largest.gov <- partydim2[(r-1),][bargpow[(r-1),]==max(bargpow[(r-1),])]      
   positions <- mapply(partyplacer,
                           p=p,
                           rule=rule[2,],
                           electionst1=elections[1,],
                           electionst2=elections[1,],
                           partydim1=partydim1[1,],
                           partydim2=partydim2[1,],
                           partydim1t2=partydim1[1,],
                           partydim2t2=partydim2[1,],                                             
                           speed=speed,
                           govdim1=govdim1[1],
                           govdim2=govdim2[1],
                           direction1=partydim1[1,]-partydim1[1,],
                           direction2=partydim2[1,]-partydim2[1,],
                           gov=gov[1,],        
                         MoreArgs=
                           list(vprefdim1=vprefdim1,
                                vprefdim2=vprefdim2,                                
                                vote=vote
                                )
                           )
 partydim1[2,] <- as.numeric(positions[1,])
 partydim2[2,] <- as.numeric(positions[2,])


#Voting 

  distances  <- t(apply(cbind(vprefdim1, vprefdim2), 1, euclidean.distancer, partydim1[2,],
                partydim2[2,]))
  vote <- apply(distances, 1, voter) #proximity voting
   #how far are voters, on avarage, from their closest party away in this party system?
  elections[2,] <- tabulate(vote, parties)
  representation[2,] <- representer(distances)

  seatsc <- t(apply(allc, 1, product, t(elections[2,]))) #seats parties contribute to coalition
  sumseats <- rowSums(seatsc) # #seats coalition has

  coalition <- rank.formater(allc, elections[2,], partydim1[2,],
                             partydim2[2,], alpha, util[1,], discount, gov[1,],
                             pprefdim1, pprefdim2, govdim1[1,], govdim2[1,],
                             govseats[1,], cur.gov=gov[2,])

  # store information
  gov[2,] <- as.matrix(allc[coalition[1],])
  govdim1[2,] <- coalition[2]
  govdim2[2,] <- coalition[3]
  govseats[2,] <- coalition[(4+(2*parties)):(3+(3*parties))]
  util[2,] <- coalition[4:(4+(parties-1))]
  bargpow[2,] <- coalition[(4+parties):((2*parties)+3)]
  caretaker[2,] <- coalition[length(coalition)]
  rule.vote.shares <- tapply(elections[2,], rule[2,], sum)/1000
  enr[2] <- eff.num.rule(rule.vote.shares)
  eccentricity[2,] <- eccentriciter(partydim1[2,], partydim2[2,],
                                        mean.vprefdim1, mean.vprefdim2)
  gov.representation[2,] <- gov.representer(vpref, govdim1[2,], govdim2[2,],
                                                sum.dist.voter.median.voter)
 


 #########
 # DRIVE #
 #########
for (r in 3:no.evo.length) {
 print(c("r=", r))
  # Party positioning for evolution

 ##  if(r>no.evo.length){ #I also disabled the computation of best rule
 ##      rule[r,] <- evolutioner(util=util[(r-memory):(r-1),],
 ##                              aspiration=aspiration,
 ##                              potential.rule=potential.rules,
 ##                              rule=rule[(r-memory):(r-1),])
 ## }
  
   ## #get position of government party with highest bargpower  
   ## dim1.largest.gov <- partydim1[(r-1),][bargpow[(r-1),]==max(bargpow[(r-1),])]
   ## dim2.largest.gov <- partydim2[(r-1),][bargpow[(r-1),]==max(bargpow[(r-1),])]      
   positions <- mapply(partyplacer,
                           p=p,
                           rule=rule[r,],
                           electionst1=elections[r-1,],
                           electionst2=elections[r-2,],
                           partydim1=partydim1[r-1,],
                           partydim2=partydim2[r-1,],
                           partydim1t2=partydim1[r-2,],
                           partydim2t2=partydim2[r-2,],                                             
                           speed=speed,
                           govdim1=govdim1[r-1],
                           govdim2=govdim2[r-1],
                           direction1=partydim1[r-1,]-partydim1[r-2,],
                           direction2=partydim2[r-1,]-partydim2[r-2,],
                           gov=gov[r-1,],        
                         MoreArgs=
                           list(vprefdim1=vprefdim1,
                                vprefdim2=vprefdim2,                                
                                vote=vote
                                )
                           )
 partydim1[r,] <- as.numeric(positions[1,])
 partydim2[r,] <- as.numeric(positions[2,])


#Voting 

  distances  <- t(apply(cbind(vprefdim1, vprefdim2), 1, euclidean.distancer, partydim1[r,],
                partydim2[r,]))
  vote <- apply(distances, 1, voter) #proximity voting
   #how far are voters, on avarage, from their closest party away in this party system?
  elections[r,] <- tabulate(vote, parties)
  representation[r,] <- representer(distances)

  seatsc <- t(apply(allc, 1, product, t(elections[r,]))) #seats parties contribute to coalition
  sumseats <- rowSums(seatsc) # #seats coalition has

  coalition <- rank.formater(allc, elections[r,], partydim1[r,],
                             partydim2[r,], alpha, util[r-1,], discount, gov[(r-1),],
                             pprefdim1, pprefdim2, govdim1[(r-1),], govdim2[(r-1),],
                             govseats[(r-1),], cur.gov=gov[r,])

  # store information
  gov[r,] <- as.matrix(allc[coalition[1],])
 if (r!=runs) {
  gov[(r+1),] <- as.matrix(allc[coalition[1],])
 }

  govdim1[r,] <- coalition[2]
  govdim2[r,] <- coalition[3]
  govseats[r,] <- coalition[(4+(2*parties)):(3+(3*parties))]
  util[r,] <- coalition[4:(4+(parties-1))]
  bargpow[r,] <- coalition[(4+parties):((2*parties)+3)]
  caretaker[r,] <- coalition[length(coalition)]
  rule.vote.shares <- tapply(elections[r,], rule[r,], sum)/1000
  enr[r] <- eff.num.rule(rule.vote.shares)
  eccentricity[r,] <- eccentriciter(partydim1[r,], partydim2[r,],
                                        mean.vprefdim1, mean.vprefdim2)
  gov.representation[r,] <- gov.representer(vpref, govdim1[r,], govdim2[r,],
                                                sum.dist.voter.median.voter)
 } #for r

for (r in (no.evo.length+1):runs) {
  # Party positioning for evolution


   rule[r,] <- evolutioner(util=util[(r-memory):(r-1),],
                              aspiration=aspiration,
                              potential.rule=potential.rules,
                              rule=rule[(r-memory):(r-1),])

  
   ## #get position of government party with highest bargpower  
   ## dim1.largest.gov <- partydim1[(r-1),][bargpow[(r-1),]==max(bargpow[(r-1),])]
   ## dim2.largest.gov <- partydim2[(r-1),][bargpow[(r-1),]==max(bargpow[(r-1),])]      
   positions <- mapply(partyplacer,
                           p=p,
                           rule=rule[r,],
                           electionst1=elections[r-1,],
                           electionst2=elections[r-2,],
                           partydim1=partydim1[r-1,],
                           partydim2=partydim2[r-1,],
                           partydim1t2=partydim1[r-2,],
                           partydim2t2=partydim2[r-2,],                                             
                           speed=speed,
                           govdim1=govdim1[r-1],
                           govdim2=govdim2[r-1],
                           direction1=partydim1[r-1,]-partydim1[r-2,],
                           direction2=partydim2[r-1,]-partydim2[r-2,],
                           gov=gov[r-1,],        
                         MoreArgs=
                           list(vprefdim1=vprefdim1,
                                vprefdim2=vprefdim2,                                
                                vote=vote
                                )
                           )
 partydim1[r,] <- as.numeric(positions[1,])
 partydim2[r,] <- as.numeric(positions[2,])


#Voting 

  distances  <- t(apply(cbind(vprefdim1, vprefdim2), 1, euclidean.distancer, partydim1[r,],
                partydim2[r,]))
  vote <- apply(distances, 1, voter) #proximity voting
   #how far are voters, on avarage, from their closest party away in this party system?
  elections[r,] <- tabulate(vote, parties)
  representation[r,] <- representer(distances)

  seatsc <- t(apply(allc, 1, product, t(elections[r,]))) #seats parties contribute to coalition
  sumseats <- rowSums(seatsc) # #seats coalition has

  coalition <- rank.formater(allc, elections[r,], partydim1[r,],
                             partydim2[r,], alpha, util[r-1,], discount, gov[(r-1),],
                             pprefdim1, pprefdim2, govdim1[(r-1),], govdim2[(r-1),],
                             govseats[(r-1),], cur.gov=gov[r,])

  # store information
  gov[r,] <- as.matrix(allc[coalition[1],])
 if (r!=runs) {
  gov[(r+1),] <- as.matrix(allc[coalition[1],])
 }

  govdim1[r,] <- coalition[2]
  govdim2[r,] <- coalition[3]
  govseats[r,] <- coalition[(4+(2*parties)):(3+(3*parties))]
  util[r,] <- coalition[4:(4+(parties-1))]
  bargpow[r,] <- coalition[(4+parties):((2*parties)+3)]
  caretaker[r,] <- coalition[length(coalition)]
  rule.vote.shares <- tapply(elections[r,], rule[r,], sum)/1000
  enr[r] <- eff.num.rule(rule.vote.shares)
  eccentricity[r,] <- eccentriciter(partydim1[r,], partydim2[r,],
                                        mean.vprefdim1, mean.vprefdim2)
  gov.representation[r,] <- gov.representer(vpref, govdim1[r,], govdim2[r,],
                                                sum.dist.voter.median.voter)
 } #for r

  ##################
  # makingdatasets #
  ##################
 
  Sprule <- matrix(NA, r*parties, 1)
  Sdim1  <- matrix(NA, r*parties, 1)
  Sdim2  <- matrix(NA, r*parties, 1)
  Spalpha <- matrix(NA, r*parties, 1)
  Spprefdim1 <- matrix(NA, r*parties, 1)
  Spprefdim2 <- matrix(NA, r*parties, 1)
  Sspeed <- matrix(NA, r*parties, 1)
  Spmemory <- matrix(memory, r*parties, 1)
  Spaspiration <- matrix(aspiration, r*parties, 1)
  Sparties <- matrix(parties, r*parties, 1)
  Srepresentation <- matrix(NA, r, 1)
  Sutil <- matrix(NA, r*parties,1)
  Sbargpow <- matrix(NA, r*parties,1)
  Sdiscount <-  matrix(discount, r*parties, 1)
  Sno.evo.length <- matrix(no.evo.length, r*parties, 1)
  Sevo.burn.in <- matrix(evo.burn.in, r*parties, 1)
  Ssd.distance <- matrix(sd.distance, r*parties, 1)
  Sseed <- matrix(run.seed, r*parties, 1)
  Svdim1.mean <- matrix(mean(as.vector(vprefdim1), na.rm=TRUE), r*parties, 1)
  Svdim2.mean <- matrix(mean(as.vector(vprefdim2), na.rm=TRUE), r*parties, 1)
  Sseats <- matrix(NA, r*parties, 1)
  Sgovmember<- matrix(NA, r*parties, 1)
  Seccentricity <- matrix(NA,r*parties, 1)
  Sgov.representation <- matrix(NA, r, 1)
  round <- rep(1:r, parties)
  party <- rep(1:parties, each=r)
  caretaker <- rep(caretaker, parties)
   
  Srepresentation <- rep(representation[1:r], parties)
  Sgov.representation <- rep(gov.representation[1:r], parties)
  Seccentricity <- rep(eccentricity[1:r], parties)
  
  for(p in 1:parties){
   Sprule[((p-1)*r+1):(p*r)] <- rbind(as.matrix(rule[(1:r),p]))
   Sdim1[((p-1)*r+1):(p*r)] <- rbind(as.matrix(partydim1[(1:r),p]))
   Sdim2[((p-1)*r+1):(p*r)] <- rbind(as.matrix(partydim2[(1:r),p]))
   Spalpha[((p-1)*r+1):(p*r)] <- alpha[p]
   Spprefdim1[((p-1)*r+1):(p*r)] <- pprefdim1[p]
   Spprefdim2[((p-1)*r+1):(p*r)] <- pprefdim2[p]
   Sspeed[((p-1)*r+1):(p*r)] <- speed[p]
   Sseats[((p-1)*r+1):(p*r)] <- rbind(as.matrix(elections[(1:r),p]))
   Sgovmember[((p-1)*r+1):(p*r)] <- rbind(as.matrix(gov[(1:r),p]))
   Seccentricity[((p-1)*r+1):(p*r)] <- rbind(as.matrix(eccentricity[(1:r),p])) 
  }
 
 dataset <- cbind(num, experiment, round, party, Sprule, Sdim1, Sdim2, Spalpha,
                   Spprefdim1, Spprefdim2, Spmemory, Spaspiration, Sparties, 
                   Sdiscount, Sno.evo.length, Sevo.burn.in, Ssd.distance, Sseed,
                   Svdim1.mean, Svdim2.mean, Sspeed, Sseats,
                   Sgovmember, Srepresentation, Seccentricity, Sgov.representation,
                   enr[1:r], caretaker)
 colnames(dataset) <- c("num", "experiment","round", "party","rule","dim1","dim2","alpha",
                        "pprefdim1", "pprefdim2", "memory", "aspiration", "parties", 
                        "discount", "no.evo.length", "evo.burn.in", "sd.distance", "seed",
                        "voter.mean.dim1", "voter.mean.dim2", "speed",
                        "size", "govmember", "representation", "eccentricity",
                        "gov.representation", "ENR", "caretaker")
 dataset <- dataset[dataset[,"round"]%in% seq((max(dataset[,"round"])-99),
                    max(dataset[,"round"])), ]
 save(dataset, file=paste0(folder, "sim_",experiment,"_", num, ".Rdata"))
} #file.exists 
}
